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We propose a new method to sum up electrostatic interactions in 2D slab geometries. It consists 
of a combination of two recently proposed methods, the 3D Ewald variant of Yeh and Berkowitz, J. 
Chem. Phys. Ill (1999) 3155, and the purely 2D method MMM2D by Arnold and Holm, to appear 
in Chem. Phys. Lett. 2002. The basic idea involves two steps. First we use a three dimensional 
summation method whose summation order is changed to sum up the interactions in a slab-wise 
„Cj ■ fashion. Second we subtract the unwanted interactions with the replicated layers analytically. The 

D ' resulting method has full control over the introduced errors. The time to evaluate the layer correction 

term scales linearly with the number of charges, so that the full method scales like an ordinary 3D 
. Ewald method, with an almost linear scaling in a mesh based implementation. In this paper we 

^s^j ' will introduce the basic ideas, derive the layer correction term and numerically verify our analytical 

results. 
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C/j ; I. INTRODUCTION 

The calculation of long range interactions due to Coulomb, gravitational, or dipolar particles is of broad interest from 
the astrophysics, the biophysics up to the solid state community. These interactions present a formidable challenge 
even to modern computers. Sophisticated methods such as fast multipole methods, tree code algorithms, Poisson grid 
solvers, or Ewald mesh methods, bring the complexity of N interacting particles down to an almost linear scaling for 
three-dimensional periodic systems. Often, however, one is interested in slab-like systems which are only periodic in 
two space dimensions and finite in the third, for example in problems involving electrolyte solutions between charged 
surfaces, proteins near charged membranes, thin films of ferrofluids, Wigner crystals, charged films, membranes, solid 
surfaces decorated with dipoles etc. 

For such systems Ewald based formulas are only slowly convergent, have mostly CHN 2 ) scalings and no "a priori" 
error estimates existEJ. East Ewald based methods have been recently put forward incl, and a non-Ewald method has 
been put forward in Rcf.ta that is based on a resummation of the force sum. However, these methods are hampered due. 
to non-controllable errors and an 0(N 2 ) scaling respectively. Recently we proposed a new method called MMM2DHlI 
which has an 0(N 5 ^ 3 ) complexity and full error control that is based on a convergence factor approach similar to 

■ MMMlI. However, this will still only allow simulations including up to a few thousand charges. There have been early 
\ attempts to use a 3D Ewald sum for these slab problems. The main idea is to fill only parts of the simulatioaJjox with 

fS| ■ charges and to leave some space empty, in an attempt to decouple the interactions in the third dimensioro'Q'Ll. Since 
each image layer is globally neutral, one hopes that their interactions decay as they become more and more distant, 
i.e. as the siae of the gap is increased. In this way one could make use of any advanced 3D Ewald implementation, 
see also Ref.E2l for a variant of this idea. 

In this paper we will follow the last suggestion and derive a term, called electrostatic layer correction (ELC), which 
subtracts the interactions due to the unwanted layers. The combination of that term with any three dimensional 
summation method with slab-wise summation order will yield the exact electrostatic energy. Since the change in the 
• summation order is done by adding a very simple term, any three dimensional summation method with the standard 
spherical summation order can be used. The new term can be evaluated easily in a time linear in the number of 
charges, hence the whole method scales like the underlying standard summation method. We develop also an error 
' formula for the maximal pairwise error in the energy and forces of the layer correction term, hence the precision of 
^ , this method can be-txped to any desired value, when used in conjunction with other error estimates for the standard 
summation metho dOE3 In the first section we will recapitulate the way how to correct the summation order via a 
modified dipole term. In the second section we will derive the layer correction term, and develop in the following 

■ section error estimates for its value. The applicability of our method will be demonstrated by a numerical analysis in 
the following section, and we end with our conclusion. 
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II. CHANGING THE SUMMATION ORDER 



We consider a system of N particles with charges qi and positions pi = (xi,yi,Zi) that reside in a box of edges 
L x L x h, where h = max^.j |z, — Zj\ is the maximal z-distance of two particles. The basic idea is to expand this 
slab system in the non-periodic z-coordinate to a system with periodicity in all three dimensions. More precisely, the 
original box of size L x L x h is placed inside a box of size L x L x L z where L z >> h sufficiently large. Then this 
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box is replicated periodically in all three dimensions. The result is a three-dimensional periodic system with empty 
space regions ("gaps") of height S := L z — h (see Fig. [I]). 5 will be called gap size in the following. 
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FIG. 1: Schematic representation of a fully periodically replicated slab system 



Since the electrostatic potential is only finite if the total system is charge neutral, the additional image layers 
(those layers above or below the original slab system) are charge neutral, too. Now let us consider the n th image 
layer which has an offset of nL z to the original layer. If nL z is large enough, each particle of charge qj at position 
(xj, yj, Zj+nL z ) and its replicas in the x, y-plane can be viewed as constituting a homogeneous charged sheet of charge 
density o~j = jj. The potential of such a charged sheet at distance z is 2it<jj\z\. Now we consider the contribution 
from a pair of image layers located at ±nL zl n > to the energy of a charge qi at position (xi, yi, Zi) in the central 
layer. Since \zj — Zi\ < nL z , we have \zj — Zi + n,L z \ = nL z + Zj — z% and \zj ~ Zi — nL z \ = nL z — Zj + Zi, and hence 
the interaction energy from those two image layers with the charge qi vanishes by charge neutrality: 

N N 

Znqi &j(\zj - Zi + nL z \ + \zj - Zi - nL z \) = \'Kq i nL z y~] (Xj = . (f) 

The only errors occurring are those coming from the approximation of assuming homogeneously charged, infinite 
sheets instead of discrete charges. This assumption should become better when increasing the distance nL z from the 
central layer. _ 

However, in a naive implementation, even large gap sizes will result in large errorsa. This is due to the order of 
summation for the three dimensional Coulomb sum, which is spherical by convention. This order implies that with 
increasing shell cutoff S the number of image shells grows faster than the number of shells of the primary layer, 



namely 0(S 3 ) versus 0(S 2 ) (see Fig. 2(a)). In other words, we include the unwanted terms faster than the actually 
wanted terms. .-Also the image layers are not really infinite charged sheets but are truncated due to the cut-off. Yeh 
and Berkowitza already suggested that this problem can be solved by changing the order of summation. Smith has 
shown that by adding to the Coulomb energy the term 

2 2irM 2 

E c = 2ttM| - — ^— , (2) 

wheije-iM = ^ qtpi is the total dipole moment, one obtains the result of a slab- wise summation instead of the spherical 
limiilia. Slab -wise summation refers to the sum X)| n |>o Ei(n), where Ei(n) denotes the energy, calculated in spherical 
summation order, resulting from the image layer with shift nL z in the z-coordinate. Technically this is the orde r 
where we first treat the original layer and then add the image layers grouped in symmetrical pairs (see Fig. |2(b)| ). 
Obviously this summation order fits much better to the charged sheet argumentation given above. Although this is 
a major change in the summation order, the difference given by Eq. (j^) is a very simple term. In fact, Smith shows 
that changes of the summation order always result in a difference that depends only on the total dipole moment. 

Applying this slab-wise summation order, Yeh and Berkowitz showed that a gap size of at least h is normally 
sufficient to obtain an moderately accurate result. Therefore the result of a standard three dimensional summation 
method plus the shape-dependent term given by Eq. (^|) , which we refer to as a slab-wise method, can be used to obtain 
a good approximation to the result for the slab geometry with the same computational effort as for the underlying 
three dimensional summation method (no matter if a simple or sophisticated method is used) . One drawback is that 
no theoretical estimates exist for the error introduced by the image layers. Therefore one might be forced to use even 
larger gaps to assure that no artifacts are produced by the image layers. One simple deducible artifact is that the 
pairwise error will be position dependant. Particles in the middle of the slab will see no effect of the image layers due 
to symmetry, and particles near the surface will encounter for the same reason the largest errors, which is definitely 
an unwanted feature for studying surface effects. Therefore averaging error measures like the commonly used RMS 
force error should not be applied without additional checks for the particles near the surfaces. 
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(a) Schematic view of the spherical summation order. S is 
the length of the box offset. 
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(b) Schematic view of the slab— wise summation order, n is 
the z offset of the box, the spherical summation order in the 
x,j/-plane is not shown. 



The other drawback is that normally the box now will have a significantly larger L z /L. But at least for Ewald 
type methods the computation time is proportional to this fraction. This is easy to see as the number of fc-space 
vectors in the z direction must be proportional to L z to maintain a fixed resolution and therefore error. It is verified 
experimentally that a gap of at least h is needed. For a cubic system h = L therefore the computation time at least 
doubles. r - in 

Nevertheless .because of the bad scaling of the known methods for slab geometries like the one by ParryoEj (0(N 2 )) 
or MMM2Dq'E| (0(iV 5 / 3 )), for particle numbers above N m 1000 using slab-wise methods is a great improvement. 



III. THE ELECTROSTATIC LAYER CORRECTION TERM 



We will now derive a term that allows to calculate the exact contribution of the image layers very efficiently, which 
we will call the electrostatic layer correction (ELC) in the following. For the following analysis there is no special 
restriction on h except for h < L z , which is true even if the L x L x L z =box is completely filled. 

The method presented here is heavily based on parts of MMM2DQ. We start with a formal definition of the 
Coulomb energy of the slab system 

oo _^ n ^ 

£= 2^ H |- ' T-j ■ ( 3 ) 

S=0 neZ 2 x {0}i,j=l ly * y 3 I 
nl+n 2 y =S 

A = diag(L, L, L z ) is a diagonal matrix describing the shape of the box. The image boxes are denoted with the 
vector n = (n x ,n y ,n z ), where n z = for now. The prime on the inner summation indicates the omission of the 
self-interaction i = j in the primary box n = (0, 0, 0) (i. e. the singular case). For the surrounding dielectric medium 
we assume vacuum boundary conditions. 

We now expand the system to a fully three-dimensional periodic system, where L z determines the period in the 
^-coordinate as in the previous section. We can rewrite the energy as 



where 



E = E S +E c + E lc , (4) 



n 2 =S 



denotes the standard three-dimensional Coulomb sum with spherical limit. To evaluate this expression one can use 
any of the efficient algorithms, starting with the classical Ewald summation up to modern methods like fast multipole 



isEl or mesh based algorithms^ 



methodsli3 or mesh based algorithmsO. E c again denotes the shape-dependent term given by Eq. (|2|) and finally 

oo TV 

^ = -\T, E E E E b . («) 

denotes the contribution of the image layers, for which we are going to derive a new expression in the following. 
We start with the expression for the energy induced by an image layer at z-offset n z ^ 0: 

oo N 

^)=4e e e lP .-T 3 +An r w 

n 2 = S 

It can be shown rigorously, although this is non-trivial, that 

Ei(n z ) = -~ lim Y Y M»£ — - . (8) 

na 2 x{ ttz } i,i=i ^* FJ 1 

This is a convergence factor approach with a convergence factor of e~^ Pi ~ Pj+An ^ . Note that this approach is exact 
only forJjaia-dimensional systems, for three-dimensional systems Eqs. (f7j) and (||) differ by a multiple of the dipolc 
moriwpalatj. 

InQcJ one can find a proof for this equation and an efficient way of calculating Ei for charge neutral systems. We 
do not want to go through the full derivation again; it consists of the application of Poisson's summation formula 
along both periodic coordinates and performing the limit j3 — > analytically. One obtains 



Ei c {n z ) = - 2 E 9iQj^>(Pi - Pj + An) , (9) 

where <f> is given by 




4 e -1-K\k\\\\z\/L 

{x,y,z)= — rj—. cos(2wk x x/L) cos(2wk y y/L) + 

e ~2-!rk y \z\/L e -27rfc x |z|/L \ 2n 

cos(2Trk y y/L) + y cos(2irk x x/L) - -^\z\ . (10) 

kx k x >Q v J 

k\\ = {k x , k y ) is a Fourier variable with integer values. <f> is an artificial pairwise potential that yields the total Coulomb 
energy and its derivative produces the pairwise forces for the periodic system. 

For now we only have a formula for the contribution of one image layer, so we still have to sum over all n z . This 
task can be performed analytically. The terms 2ir\z\/L 2 can be omitted since they are exactly the homogeneous sheet 
potential and we have seen before that this cancels out for charge neutral systems (see Eq. (|l|)). 

The summation over n z of the remaining sums over p and q is fairly easy to perform using the geometric series 
(as these sums are absolutely convergent, exchanging the summation over n z and the summations over (k X: ky^) is 
possible). Combining the terms for ±n z again we obtain 

Eic = y qiqjip(Pi - Pj + An) , (11) 



where 



4 ^-^ cosh(2irk\\Zij / L) 

4>{x, y, z) =- y fc ^fcuWL _ ^ cos^irkxXij/L) cas{2-Kk y y ij / L) + 

2 ^ coah(2irk x Zij / L) cos(2irk x Xij / L) 

L ^ k x (e 2 * k * L '/ L - 1) + ^ 

2 ^-^ cosh(27tky Zij / L) cos(2Trk y yij/ L) 

L ^ L(e M »^/ L - 1) ' 

k v >0 yy ' 
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The forces can be obtained from that by simple differentiation since the sums are absolutely convergent. Although 
the form in Eq.(|l^) has a much better convergence than the original form in Eq.(Q), its main advantage is a linear 
computation time with respect to the number of particles N. To see this, the equation has to be rewritten using the 
addition theorems for the cosine and the hyperbolic cosine. For each fcy one first calculates the sixteen terms 



X(c/s,c/s,c/s) = cosh/ sinh(27rfc||Z i /L)cos/ sin(2nk x Xi/L) cos / sm(2wk y yi/L) , 
»=i 

N 

X(x,c/s,c/s) = X! * cosh / sili H2nk x Zi/L) cos / sin(2irk x x l / L) , (13) 

i—l 
N 

X(y,c/a,c/s) = y^gi cosh/ sinh(2nky Zi / L) cos/ sm(2nk y y l /L) , 



where the indices in the obvious way determine which of the functions cosine (hyperbolicus) or sinus (hyperbolicus) 
are used. Then we evaluate 

Elc = T (p 2^L z /L _ m (X(ccc) + ^(csc) + Xfccs) + xfcss) ~ 

k m ,k y >Q ^ ' II 

X(scc) X(ssc) X(scs) X(sss)^j 

2 i • < 14 > 



1 ( 1 i 2 _ 2 _2 \ , 

^ e 2Trk x L z /L _ ]^£ x V^( ICC ) X(xcs) X(a; S c) X( XS s) J "+ 



fei>0 

2 



1 / 2 , 2 _2 _ 2 \ 

( e 2irk y L z /L _ ^%cc) + Afj/cs) X(j,sc) X(yss) J 



L ^ ( e ^k y L z /L _ i\ k 
k y >0 y ' 

Similar expansions using the same sixteen terms can also be found for the forces. Obviously this has linear 
computation time with respect to the number of particles, as the only summations over all the particles occur in 
the x*- But up to now there is still the infinite summation over ku. So the next task is to derive an estimate for the 
error induced by the replacement of the infinite sum by a finite one. 

IV. ERROR ESTIMATES 

Since Ei c is written as sum over an alternative potential ?/>, it is reasonable to derive an upper bound for the 
error from the calculation of ij) only with a finite cutoff. From this upper bound, crude estimates for other error 
measures such as the RMS (root-mean-square) force error can be derived. Again these error estimates are taken from 
MMM2Ena. As we will show later, the error distribution is not uniform along the z-axis. The error gets maximal 
for particles near the borders. Since these particles will normally be those of special interest, the maximal pairwise 
error should be some magnitudes smaller than the thermal noise. 

While the error bounds for MMM2D were only used to tunc the algorithm, the error estimates for Ei c can also 
be used to obtain an error bound for the slab-wise method from Ref.EI, and hence one can determine "a priori" the 
necessary gap size to reach a preset precision. Therefore we also have to deal with small cutoffs, especially the case 
when no terms of Ei c are added. 

First we choose a cutoff R > 1 and then evaluate Ei c only over the area 

T R = {{k X , ky) S I? I k X , ky > 0, (k X - l) 2 + (ky - 1 ) 2 < R 2 } U 

{(k x ,0)eZx{0}\k x <R} (15) 

{(0,ky) G {0} X 1\ky < R} . 
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The three sets correspond to the three sums in Eq.(|12|). Therefore we actually evaluate 

4 v-^ cosh(2n k\\ Zij / L) 

Eic=-j- 2^ fc|| ( e 2*-fcuL,/L _ ^ cos{2irk x Xij/L) cos(27rfc y ^ 3 /L) 



1 ^ k x (e 2 " k * L >/ L - 1) + 

0<k x <R xy 1 

2 x ^ cosh(27rk y Zij / L) cos(2Trk y yij / L) 



E 



L ^ kJe 2vk y L ^ L - 1) 

0<ky<R Vy ' 

Tr may look more complicated then necessary. But this form enables us to find a rigorous upper bound for the 
error. An upper bound for the absolute value of the summands is 



cash(2irk\\Zij/L) 

cos(Z7Tk x Xij j L) cos(2irk y yij / L) 



< 



-2*k„Lz/L cosh(2^fc||Z u -/L) < c _ 2nknLz/L cosh(27rfc||fe/L) 



Of course because the cosine hyperbolicus is monotonous, one could use any larger value for h. This is for example 
necessary in a priori estimations. Using this we find the upper bound for the maximal pairwise error in the potential 
by a simple approximation of the sums by integrals as 

_ 1/2 + (ttR)- 1 /cxp(2-KRh/L) exp(-2wRh/L)\ 

TE '— e 2irRL z /L _ I Lz-h + Lz + h J ' 

Details can again be found By an analogous derivation one can also find an upper bound on the maximal 

pairwise error on each of the three force components as 



TF 




exp(2irRh/L) 
(L z -h) + 

\ (18) 
exp(-2nRh/L) \ 



(L z + h) ) 



This is also an weaker bound for the potential. In other words, the maximal pairwise error on the forces is larger 
than the error in the potential. For R = 1 one obtains an overall estimate of the magnitude of the contribution of the 
image layers, i. e. an error estimate for the slab-wise methods. 

Note that Eq.(|lj]) shows that the error in the potential or the force for a single particle will be largest if it is located 
near the gap, since there \zij\ will be maximal. This effect will increase with increasing R. Therefore when using the 
layer correction one must apply non-averaging error estimates such as our maximal pairwise error. Averaging error 
estimates such as the RMS force error might be misleading about the error on the particles of interest. 

Moreover the error will decrease exponentially with the distance from the gap. Since the particles near the gap 
(i.e. the surface) are normally of special interest in simulations with slab geometry, averaging error measures like the 
RMS force error might be misleading about the effect of these errors. 

All our error estimates show that the error drops exponential both with R and L z /L. The decay in R means that 
it is easy to achieve high accuracies with our layer correction formula, while the decay in L z /L shows that slab- wise 
methods can achieve good accuracies without increasing L z /L too much. 

Although we do not encourage using the RMS error measure as we explained above, we still want to give an upper 
bound on the average error in Ei c and the RMS force error. We assume that the pairwise potentials of the different 
particle pairs are independent identically distributed random variables, which is true for homogenous random systems 
and normally a good assumption otherwise. Of course the self interaction of the particles, i. e. qftjj(Q, 0, 0) has to be 
omitted. Let erg be the variance of this random variable, then it is easy to see that oe < t e . Using this one can show 
that 

(e 1c -J2 9i V»(0, 0, 0) \ < Q 2 ^ < Q 2 t e , (19) 
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where Q := J2i Qi- 

Similarly we define ap as the variance of the forces measured in the Euclidian norm. Then because of the component- 
wise maximal error estimate for the force, we have ap < 3rj=>, and one obtains 

(^Jl/NjpA^\j < V3Q 2 /VNt f , (20) 

where AF/ C denotes the error in the layer correction force on particle i. Note that both estimates (|l9|) and (|2(]) are 
much larger than the real error as one expects a to be much smaller the maximal error (about 2 — 4 magnitudes). 



V. NUMERICAL DEMONSTRATION 
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FIG. 2: The measured error in the potential and the force of the ELC-term versus the estimates te, tf for different cutoff 
radii R 



In this section we want to show results from our implementation of the layer correction term (ELC). First we want 
to show that our maximal pairwise error bounds are correct. To this aim we place two particles randomly in a box of 
size 1 x 1 x 0.8, so that we leave a gap of 5 = 0.2 in a box of dimensions lxlxl. Fig. || shows the maximal potential 
and force error that occurred during 10000 evaluations and our estimates te and Tp. Moreover we have included the 
result for a particle pair with a relative position of (0,0, ft.), the worst case position. For such a position the error 
estimate is exact up to the approximation of the sum by an integral. As "exact" force we used R = 30. 

One can see that the maximal error estimates are always above the measured deviations. Even after 10000 random 
evaluations the maximum error is considerably lower than for the special pair with relative position (0,0, h), and 
the error is still not very smooth. This effects are due to the Fourier representation with exponentially decaying 
coefficients, which makes the worst case error extremely rare. But in a real simulation the particle distribution is not 
necessarily homogeneous and to be on the safe side one has to deal also with the rare worst case error. Nevertheless 
Fig. H shows that the error coming from the image layers can be strictly controlled. 
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(a) Computation times for the ELC term for different system sizes 
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Next we investigate the comp utati on times of our implementation of the ELC term. For different gap sizes 5 = 
0.05,0.1 and 0.2 we show in Fig. 3(a) the computation time Tcpu for different numbers of particles N. The systems 
were again consisting of uniformly randomly distributed particles of charges ±1. The maximal pairwise error was 
fixed to be 10~ 4 . The times are averages over 10 runs on a Compaq XP1000. The computation time for the same 
system consisting of 1000 charges using P 3 M is » 330ms for a typical RMS force error of 10~ 4 . Therefore even the 
small gap of 0.05 gives just the same order of computat ion t ime. For the normally used gap sizes of « 0.2L Z the 
computation time is negligible compared to P 3 M. In Fig. |3(b)| we show the RMS force errors that occurred. One can 
see the predicted Q 2 /y/N behavior. We also show the theoretical upper bound V3Q 2 /\/Nt f for R = 13 and 8 = 0.2, 
which is considerably above as expected. 
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FIG. 4: Error distribution of the layer correction along the z-axis for 100 random systems with 100 particles 



As the last important fact we demonstrate that the errors for the layer correction indeed are maximal near the gap 
(i.e. near the surface). For Fig. ^ 100 particles were put 100 times randomly in a box with a gap of 6 = 0.1. Then 
for every particle the magnitude of the layer correction for R = 40, which is a good approximation to the full Ei c , 
and the difference in the layer correction between R = 5 and R = 40 was drawn against the z-coordinate. Clearly 
the error always is largest near the gap. This effect increases with increasing R, which is easy to understand from 
the error formula. Therefore the full RMS error of the system might be completely misleading about the effect the 
errors have on the particles near the gap as we mentioned before. Nevertheless the figure shows that Ei c with R = 5 
reduces the error near the surfaces by a factor of ~ 100. 



VI. CONCLUSION 



We have derived a term called ELC to efficiently calculate the contribution of the image layers in three dimensionally 
periodically replicated slab systems. ELC scales as the number N of particles and has a rigorous error bound. 
Moreover this error bound can be used to estimate the size of the image layer contribution and therefore gives a 
bound on the error introduced by slab-wise methods as proposed by Yeh and Berkowitz. We have found that the 
error for these methods decays exponentially in L z /L. However, the errors are not uniformly distributed over the 
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slab, namely they are worst at the surfaces of the slabs. This strongly suggests to restrict the maximal pairwise error 
instead of the usually ass-tuned RMS-errors. 

In a forthcoming paperE^I we will focus on the application of ELC to the standard Ewald method and to P 3 M. We 
will show how the error formulas of Kolafa and PerramEj have to be adapted to allow non-cubic simulation boxes 
which is essential for using ELC with R = 1, i. e. a slab-wise method. For all combinations we present numerical 
results which allow easily to decide which method is optimal for use in a real simulation. 
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